library(foreign)
x <- read.dta("http://www-personal.umich.edu/~jowei/fema/Chen.Individual.Data.dta")
install.packages("xtable")
require(xtable)

## Subset the main dataset to desired groups

## Female
x.female <- subset(x, x$male==0)

## Male
x.male <- subset(x, x$male==1)

## African Americans
x.black <- subset(x, x$black==1)

## Non-African Americans
x.nonblack <- subset(x, x$black==0)

## Table 1, Model 1: Treatment Effect on Male Democrats 
male.democrats <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, 
                              family=binomial(link='logit'), 
                              data=x.male[x.male$Party=="DEM" & x.male$date<=41104,]))
xtable(male.democrats)

## Table 1, Model 2: Treatment Effect on Male Republicans
male.republcians <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, 
                              family=binomial(link='logit'), 
                              data=x.male[x.male$Party=="REP" & x.male$date<=41104,]))
xtable(male.republicans)

## Table 1, Model 4: Treatment Effect on Male Democrats with Covariates 
male.democrats.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+medhomevalue10k+
              (Gender=="M")+black+medincome10k+as.factor(CountyCode), family=binomial(link='logit'), 
            data=x.male[x.male$Party=="DEM" & x.male$date<=41104,]))
xtable(male.democrats.cov)

## Table 1, Model 5: Treatment Effect on Male Republicans with Covariates 
male.republicans.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+
              medhomevalue10k+(Gender=="M")+black+medincome10k+as.factor(CountyCode), 
            family=binomial(link='logit'), data=x.male[x.male$Party=="REP" & x.male$date<=41104,]))
xtable(male.republicans.cov)

## Table 1, Model 1: Treatment Effect on Female Democrats
female.democrats <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, 
                              family=binomial(link='logit'), 
                              data=x.female[x.female$Party=="DEM" & x.female$date<=41104,]))

## Table 1, Model 2: Treatment Effect on Female Republicans
female.republicans <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, 
                                family=binomial(link='logit'), 
                                data=x.female[x.female$Party=="REP" & x.female$date<=41104,]))
xtable(female.republicans)

## Table 1, Model 4: Treatment Effect on Female Democrats with Covariates
female.democrats.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+medhomevalue10k+
              (Gender=="M")+black+medincome10k+as.factor(CountyCode), family=binomial(link='logit'), 
            data=x.female[x.female$Party=="DEM" & x.female$date<=41104,]))
xtable(female.democrats.cov)

## Table 1, Model 5: Treatment Effect on Female Republicans with Covariates
female.republicans.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+
              medhomevalue10k+(Gender=="M")+black+medincome10k+as.factor(CountyCode), 
            family=binomial(link='logit'), data=x.female[x.female$Party=="REP" & x.female$date<=41104,]))
xtable(female.republicans.cov)

## Table 1, Model 1: Treatment Effect on African American Democrats
black.democrats <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, family=binomial(link='logit'), 
            data=x.black[x.black$Party=="DEM" & x.black$date<=41104,]))
xtable(black.democrats)

## Table 1, Model 2: Treatment Effect on African American Republicans
black.republicans <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, family=binomial(link='logit'), 
                               data=x.black[x.black$Party=="REP" & x.black$date<=41104,]))
xtable(black.republicans)

## Table 1, Model 4: Treatment Effect on African American Democrats with Covariates
black.democrats.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+
              medhomevalue10k+(Gender=="M")+black+medincome10k+as.factor(CountyCode), 
            family=binomial(link='logit'), data=x.black[x.black$Party=="DEM" & x.black$date<=41104,]))
xtable(black.democrats.cov)

## Table 1, Model 5: Treatment Effect on African American Republicans with Covariates 
black.republicans.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+
              medhomevalue10k+(Gender=="M")+black+medincome10k+as.factor(CountyCode), 
            family=binomial(link='logit'), data=x.black[x.black$Party=="REP" & x.black$date<=41104,]))
xtable(black.republicans.cov)

## Table 1, Model 1: Treatment Effect on Non African American Democrats
nonblack.democrats <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, family=binomial(link='logit'), 
            data=x.nonblack[x.nonblack$Party=="DEM" & x.nonblack$date<=41104,]))
xtable(nonblack.democrats)

## Table 1, Model 2: Treatment Effect on Non African American Republicans
nonblack.republicans <- summary(glm(Turnout2004~Turnout2002+Aid_awarded, family=binomial(link='logit'), 
                                  data=x.nonblack[x.nonblack$Party=="REP" & x.nonblack$date<=41104,]))
xtable(nonblack.republicans)

## Table 1, Model 4: Treatment Effect on Non African American Democrats with Covariates
nonblack.democrats.cov <- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+
              medhomevalue10k+(Gender=="M")+black+medincome10k+as.factor(CountyCode), 
            family=binomial(link='logit'), data=x.nonblack[x.nonblack$Party=="DEM" & x.nonblack$date<=41104,]))
xtable(nonblack.democrats.cov)

## Table 1, Model 5: Treatment Effect on Non African American Republicans with Covariates 
nonblack.republicans.cov<- summary(glm(Turnout2004~Turnout2002+Aid_awarded+mph+mphSQUARED+age+ageSQUARED+
              medhomevalue10k+(Gender=="M")+black+medincome10k+as.factor(CountyCode), 
            family=binomial(link='logit'), data=x.nonblack[x.nonblack$Party=="REP" & x.nonblack$date<=41104,]))
xtable(nonblack.republicans.cov)

## Simulation 

install.packages("Zelig")
require(Zelig)

## Male Republicans 

zelig.out.rep.male <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                   data=x.male[x.male$Party=="REP" & x.male$date<=41104,], model="logit")
summary(zelig.out.rep.male)

x.evs0.rep.male <- c(1,mean(x.male[x.male$Party=="REP" & x.male$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.rep.male) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.rep.male <- as.data.frame(t(x.evs0.rep.male))
x.evs1.rep.male <- c(1,mean(x.male[x.male$Party=="REP" & x.male$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.rep.male) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.rep.male <- as.data.frame(t(x.evs1.rep.male))

zelig.sim0.rep.male <- sim(zelig.out.rep.male, x=x.evs0.rep.male, x1=x.evs1.rep.male, num=10000)
summary(zelig.sim0.rep.male)
plot(zelig.sim0.rep.male)

## Male Democrats 

zelig.out.dem.male <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                            data=x.male[x.male$Party=="DEM" & x.male$date<=41104,], model="logit")
summary(zelig.out.dem.male)

x.evs0.dem.male <- c(1,mean(x.male[x.male$Party=="DEM" & x.male$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.dem.male) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.dem.male <- as.data.frame(t(x.evs0.dem.male))
x.evs1.dem.male <- c(1,mean(x.male[x.male$Party=="DEM" & x.male$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.dem.male) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.dem.male <- as.data.frame(t(x.evs1.dem.male))

zelig.sim0.dem.male <- sim(zelig.out.dem.male, x=x.evs0.dem.male, x1=x.evs1.dem.male, num=10000)
summary(zelig.sim0.dem.male)
plot(zelig.sim0.dem.male)

## Female Republicans 

zelig.out.rep.female <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                            data=x.female[x.female$Party=="REP" & x.female$date<=41104,], model="logit")
summary(zelig.out.rep.female)

x.evs0.rep.female <- c(1,mean(x.female[x.female$Party=="REP" & x.female$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.rep.female) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.rep.female <- as.data.frame(t(x.evs0.rep.female))
x.evs1.rep.female <- c(1,mean(x.female[x.female$Party=="REP" & x.female$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.rep.female) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.rep.female <- as.data.frame(t(x.evs1.rep.female))

zelig.sim0.rep.female <- sim(zelig.out.rep.female, x=x.evs0.rep.female, x1=x.evs1.rep.female, num=10000)
summary(zelig.sim0.rep.female)
plot(zelig.sim0.rep.female)

## Female Democrats 

zelig.out.dem.female <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                              data=x.female[x.female$Party=="DEM" & x.female$date<=41104,], model="logit")
summary(zelig.out.dem.female)

x.evs0.dem.female <- c(1,mean(x.female[x.female$Party=="DEM" & x.female$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.dem.female) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.dem.female <- as.data.frame(t(x.evs0.dem.female))
x.evs1.dem.female <- c(1,mean(x.female[x.female$Party=="DEM" & x.female$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.dem.female) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.dem.female <- as.data.frame(t(x.evs1.dem.female))

zelig.sim0.dem.female <- sim(zelig.out.dem.female, x=x.evs0.dem.female, x1=x.evs1.dem.female, num=10000)
summary(zelig.sim0.dem.female)
plot(zelig.sim0.dem.female)

## African American Republicans 

zelig.out.rep.black <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                              data=x.black[x.black$Party=="REP" & x.black$date<=41104,], model="logit")
summary(zelig.out.rep.black)

x.evs0.rep.black <- c(1,mean(x.black[x.black$Party=="REP" & x.black$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.rep.black) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.rep.black <- as.data.frame(t(x.evs0.rep.black))
x.evs1.rep.black <- c(1,mean(x.black[x.black$Party=="REP" & x.black$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.rep.black) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.rep.black <- as.data.frame(t(x.evs1.rep.black))

zelig.sim0.rep.black <- sim(zelig.out.rep.black, x=x.evs0.rep.black, x1=x.evs1.rep.black, num=10000)
summary(zelig.sim0.rep.black)
plot(zelig.sim0.rep.black)

## African American Democrats

zelig.out.dem.black <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                             data=x.black[x.black$Party=="DEM" & x.black$date<=41104,], model="logit")
summary(zelig.out.dem.black)

x.evs0.dem.black <- c(1,mean(x.black[x.black$Party=="DEM" & x.black$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.dem.black) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.dem.black <- as.data.frame(t(x.evs0.dem.black))
x.evs1.dem.black <- c(1,mean(x.black[x.black$Party=="DEM" & x.black$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.dem.black) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.dem.black <- as.data.frame(t(x.evs1.dem.black))

zelig.sim0.dem.black <- sim(zelig.out.dem.black, x=x.evs0.dem.black, x1=x.evs1.dem.black, num=10000)
summary(zelig.sim0.dem.black)
plot(zelig.sim0.dem.black)

## Non African American Republicans 

zelig.out.rep.nonblack <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                             data=x.nonblack[x.nonblack$Party=="REP" & x.nonblack$date<=41104,], model="logit")
summary(zelig.out.rep.nonblack)

x.evs0.rep.nonblack <- c(1,mean(x.nonblack[x.nonblack$Party=="REP" & x.nonblack$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.rep.nonblack) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.rep.nonblack <- as.data.frame(t(x.evs0.rep.nonblack))
x.evs1.rep.nonblack <- c(1,mean(x.nonblack[x.nonblack$Party=="REP" & x.nonblack$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.rep.nonblack) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.rep.nonblack <- as.data.frame(t(x.evs1.rep.nonblack))

zelig.sim0.rep.nonblack <- sim(zelig.out.rep.nonblack, x=x.evs0.rep.nonblack, x1=x.evs1.rep.nonblack, num=10000)
summary(zelig.sim0.rep.nonblack)
plot(zelig.sim0.rep.nonblack)

## Non African American Democrats

zelig.out.dem.nonblack <- zelig(Turnout2004~Turnout2002+Aid_awarded, 
                             data=x.nonblack[x.nonblack$Party=="DEM" & x.nonblack$date<=41104,], model="logit")
summary(zelig.out.dem.nonblack)

x.evs0.dem.nonblack <- c(1,mean(x.nonblack[x.nonblack$Party=="DEM" & x.nonblack$date<=41104,]$Turnout2002, na.rm=T),0)
names(x.evs0.dem.nonblack) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs0.dem.nonblack <- as.data.frame(t(x.evs0.dem.nonblack))
x.evs1.dem.nonblack <- c(1,mean(x.nonblack[x.nonblack$Party=="DEM" & x.nonblack$date<=41104,]$Turnout2002, na.rm=T),1)
names(x.evs1.dem.nonblack) <- c("Intercept", "Turnout2002", "Aid_awarded")
x.evs1.dem.nonblack <- as.data.frame(t(x.evs1.dem.nonblack))

zelig.sim0.dem.nonblack <- sim(zelig.out.dem.nonblack, x=x.evs0.dem.nonblack, x1=x.evs1.dem.nonblack, num=10000)
summary(zelig.sim0.dem.nonblack)
plot(zelig.sim0.dem.nonblack)

